#
#This file is intended to calibrate by hand
#
########### 
#PREAMBLE ====
###############
library(SQUAREM)
library(lfe)
library(fixest)
library(doParallel)
library(doRNG)
library(tictoc)
library(data.table)
library(matrixStats)
library(HeadR)
#if you opened the session by doubleclicking the R project in the BLP folder then switch directories.
setwd("~/Dropbox/BLP_REStat/Replication/Rfiles")
rm(list=ls())
source("CF_MMX_functions.R")
source("CF_MMX_oly_functions.R")
source("CF_MMX_MLOG_functions.R")
source("CF_MMX_richsubs_functions.R")
source("CF_MMX_calibration_functions.R")
library(Rcpp)
Rcpp::sourceCpp("../Cppfiles/crossDelta.cpp")
#
#Parallel processing setup
cl <- detectCores()-1
registerDoParallel(cores=cl)
getDoParWorkers()
seedn <- 5
MAXIT <- 500 # our standard max for FPI in the eqbm() function
#
dol.units <- 1000 # count in 000s of dollars (for income variables)
#
#############################
#Parameters common to all====
#############################
seedn <- 777 # 666
n.hh <- 1000 # households
ntb <- 0 # this is the ntb part of tau, not needed anymore 
tar <- 0.1 # this is the tariff: 10%
dist.ave.pars <- c(0.1,0.1)
n.countries <- 3
n.models <- 90
models.perfirm <- 10
n.x <- 4  # dimensions of quality (number of measured model attributes) 
U0 <- 1  # outside good, in this case there will be constant x0 =1, 
#U0 <- 0 # no outside good, x1..x{n.x}
OGpct <- 90
GINI <- FALSE
sigma.ly.factor <- 1
n.reps <- 1000
nv <- 3 #nrow(settings.var)
# settings based on data that don't change
source("CF_MMX_settings_fixed.R")
# targeted settings
source(paste0("CF_MMX_settings_calibration_","OG",OGpct,".R"))
#
############################
#Calibration starts here====
############################
#
###
#### OG90 ====
###
#
# MCES ====
#
#Panel B: MC approx of  Multiproduct Oly 
mix.var <- "MCES/MC"
n.firms <- n.models/models.perfirm
pan.var <- "b"
#Settings: some fixed (in settings file), and the others calibrated#
settings.var <- cbind(b.settings.calib,fixed.settings.mces)
saveRDS(settings.var,file=paste0("../Tables/CF_MM/",mix.var,"/CF_settings_",pan.var,"OG_",OGpct,".rds"))
#
#Implement 
#generate needed data
co.own.mat <- make.ownership(n.firms,n.models)
set.seed(seedn)
monte.func.mm.calib(repi=1,settings=settings.var,v=3,market=1,cf.meth="EHA")
#simulate
tic()
monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
              monte.rep.mm.calib(parallel=T,settings=settings.var,v,market=1)})
toc()
#moments and output
monte.moments <- monte.moments.mm(monte.allv,round.moments = TRUE)
monte.moments
monte.moments.output(monte.moments.data = monte.moments)
#
rm(monte.allv,monte.moments,co.own.mat,settings.var)
#
# MLOG 
#
#Panel A: MC approx of MC DGP 
mix.var <- "MLOG/MC"
n.firms <- n.models
pan.var <- "a"
#Settings: some fixed (in settings file), and the others calibrated#
settings.var <- cbind(a.settings.mlog.calib,fixed.settings.mlog)
saveRDS(settings.var,file=paste0("../Tables/CF_MM/",mix.var,"/CF_settings_",pan.var,"OG_",OGpct,".rds"))
#Implement 
#generate non-random data
co.own.mat <- make.ownership(n.firms,n.models)
set.seed(seedn)
#eqbm.shares <- eqbm.calib(settings=settings.var,v=3,market=1)
monte.func.mm.calib(repi=1,settings=settings.var,v=3,market=1,cf.meth="EHA")
#simulate
tic()
monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
  monte.rep.mm.calib(parallel=T,settings=settings.var,v,market=1)})
toc()
#moments and output
monte.moments <- monte.moments.mm(monte.allv,round.moments = TRUE)
monte.moments
monte.moments.output(monte.moments.data = monte.moments)
#
rm(monte.allv,monte.moments,co.own.mat,settings.var)
#
#Panel B: MC approx of OLY DGP 
mix.var <- "MLOG/MC"
n.firms <- n.models/models.perfirm
pan.var <- "b"
#Settings: some fixed (in settings file), and the others calibrated#
settings.var <- cbind(b.settings.mlog.calib,fixed.settings.mlog)
saveRDS(settings.var,file=paste0("../Tables/CF_MM/",mix.var,"/CF_settings_",pan.var,"OG_",OGpct,".rds"))
#Implement 
#generate non-random data
co.own.mat <- make.ownership(n.firms,n.models)
set.seed(seedn)
monte.func.mm.calib(repi=1,settings=settings.var,v=1,market=1,cf.meth="EHA")
#simulate
tic()
monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
  monte.rep.mm.calib(parallel=T,settings=settings.var,v,market=1)})
toc()
#moments and output
monte.moments <- monte.moments.mm(monte.allv,round.moments = TRUE)
monte.moments
monte.moments.output(monte.moments.data = monte.moments)
#
rm(monte.allv,monte.moments,co.own.mat,settings.var)
#


 ###
 #### OG10 ====
 ###
 #reset the settings
 OGpct <- 10
 source(paste0("CF_MMX_settings_calibration_","OG",OGpct,".R"))
 #
 #MCES
 #
 #Panel B: MC approx of  Multiproduct Oly 
 mix.var <- "MCES/MC"
 n.firms <- n.models/models.perfirm
 pan.var <- "b"
 #Settings: some fixed (in settings file), and the others calibrated#
 settings.var <- cbind(b.settings.calib,fixed.settings.mces)
 saveRDS(settings.var,file=paste0("../Tables/CF_MM/",mix.var,"/CF_settings_",pan.var,"OG_",OGpct,".rds"))
 #
 #Implement 
 #generate needed data
 co.own.mat <- make.ownership(n.firms,n.models)
 set.seed(seedn)
 monte.func.mm.calib(repi=1,settings=settings.var,v=3,market=1,cf.meth="EHA")
 #simulate
 tic()
 monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
   monte.rep.mm.calib(parallel=T,settings=settings.var,v,market=1)})
 toc()
 #moments and output
 monte.moments <- monte.moments.mm(monte.allv,round.moments = TRUE)
 monte.moments
 monte.moments.output(monte.moments.data = monte.moments)
 #
 rm(monte.allv,monte.moments,co.own.mat,settings.var)
 #
 #Panel C: OLY approx of  Multiproduct Oly : not needed anymore (except settings saving) since the target is now price elas and not eta.
 mix.var <- "MCES/OLY"
 n.firms <- n.models/models.perfirm
 pan.var <- "c"
 # #Settings: some fixed (in settings file), and the others calibrated#
 settings.var <- cbind(c.settings.calib,fixed.settings.mces)
 saveRDS(settings.var,file=paste0("../Tables/CF_MM/",mix.var,"/CF_settings_",pan.var,"OG_",OGpct,".rds"))
 #
 #MLOG
 #
 #Panel B of: MC approx of OLY DGP 
 mix.var <- "MLOG/MC"
 n.firms <- n.models/models.perfirm
 pan.var <- "b"
 #Settings: some fixed (in settings file), and the others calibrated#
 settings.var <- cbind(b.settings.mlog.calib,fixed.settings.mlog)
 saveRDS(settings.var,file=paste0("../Tables/CF_MM/",mix.var,"/CF_settings_",pan.var,"OG_",OGpct,".rds"))
 #Implement 
 #generate non-random data
 co.own.mat <- make.ownership(n.firms,n.models)
 set.seed(seedn)
 monte.func.mm.calib(repi=1,settings=settings.var,v=1,market=1,cf.meth="EHA")
 #simulate
 tic()
 monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
   monte.rep.mm.calib(parallel=T,settings=settings.var,v,market=1)})
 toc()
 #moments and output
 monte.moments <- monte.moments.mm(monte.allv,round.moments = TRUE)
 monte.moments
 monte.moments.output(monte.moments.data = monte.moments)
 #
 rm(monte.allv,monte.moments,co.own.mat,settings.var)
 
 ###
 #### OG50 ====
 ###
 #reset the settings
 OGpct <- 50
 source(paste0("CF_MMX_settings_calibration_","OG",OGpct,".R"))
 #
 #MLOG
 #
 #Panel B of: MC approx of OLY DGP 
 mix.var <- "MLOG/MC"
 n.firms <- n.models/models.perfirm
 pan.var <- "b"
 #Settings: some fixed (in settings file), and the others calibrated#
 settings.var <- cbind(b.settings.mlog.calib,fixed.settings.mlog)
 saveRDS(settings.var,file=paste0("../Tables/CF_MM/",mix.var,"/CF_settings_",pan.var,"OG_",OGpct,".rds"))
 #Implement 
 #generate non-random data
 co.own.mat <- make.ownership(n.firms,n.models)
 set.seed(seedn)
 monte.func.mm.calib(repi=1,settings=settings.var,v=1,market=1,cf.meth="EHA")
 #simulate
 tic()
 monte.allv <- as.data.table(foreach(v=1:nv, .combine = 'rbind')  %do% {cat(paste(mix.var,pan.var,": Starting setting ",v,"\n"))
   monte.rep.mm.calib(parallel=T,settings=settings.var,v,market=1)})
 toc()
 #moments and output
 monte.moments <- monte.moments.mm(monte.allv,round.moments = TRUE)
 monte.moments
 monte.moments.output(monte.moments.data = monte.moments)
 #
 rm(monte.allv,monte.moments,co.own.mat,settings.var)
 
 


